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FAP 

ENTRY  DPA 
DPA  CLA#  1,4 
FAD#  3 *4 
STO  TEMP 
XCA 

UFA#  2»4 
UFA#  4,4 
FAD  TEMP 
STO#  5.4 
STO#  6.4 
TRA  7.4 
TEMP  HTR  0 
END 


FAP 

ENTRY  DPS 
DPS  CLA#  1.4 
FSB#  3.4 
STO  TEMP 
XCA 

UFA#  2.4 
UFS#  4.4 
FAD  TEMP 
STO#  5.4 
STO#  6.4 
TRA  7.4 
TEMP  HTR  0 
END 
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FNOL2,  A  FORTRAN  (IBM7090)  SUBROUTINE 
FOR  THE  SOLUTION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS 
WITH  AUTOMATIC  ADJUSTMENT  OF  THE  INTERVAL  OF  INTEGRATION 

1.  INTRODUCTION 

A.  Background 

Many  techniques  are  available  for  the  approximate  solution  of 
ordinary  differential  equations  by  nvunerical  methods.  The  subroutine  des¬ 
cribed  in  this  report  exploits  the  4th  order  Runge-Kutta  and/or  Adams-Moulton 
methods  which  have  been  quite  popular  at  NOL. 

Several  years  ago  (1959),  NOLI,  reference  (2),  one  of  the  predecessors 
of  the  present  subroutine, was  introduced.  It  was  used  with  great  success  at 
NOL  and  the  techniques  it  used  for  automatic  linkage,  termination  and  output 
were  carried  over  in  FNOLl,  reference  (1),  the  forerunner  of  FN0L2. 

FNOLl  was  written  to  overcome  the  following  two  limitations  of 
reference  (2)&  1.  the  inability  to  change  step  size  efficiently,  and 

2.  the  difficulty  of  modifying  NOLI  for  special  purposes.  Consequently, 

FNOLl  was  written  in  FORTRAN  with  optional  automatic  adjustment  of  step  size. 
Unfortunately,  the  double  precision  summation  technique  of  reference  (2)  was 
not  retained. 

B.  Description  of  the  Package 
1.  Special  Features 

a.  The  addition  of  double  precision  arithmetic  in  key  locations 
has  made  the  accuracy  of  FN0L2  comparable  to  NOLI,  reference 
(2). 
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b.  The  user  may  specify  whether  the  automatic  adjustment  of  the 
integration  interval  (step  size)  is  to  be  a  fvinction  of 
RELATIVE  or  ABSOLUTE  TRUNCATION  ERROR.* 

2 .  Coding  Information 

a.  FN0L2  is  a  Fortran  Program  which  requires  4004^  cells  of 
memory . 

b.  No  COMMON  storage  is  used. 

C.  Objectives  of  the  Report 

1.  To  provide  the  reader  with  a  complete  description  of  the  use 
of  FN0L2.  See  Sections  III  and  IV. 

2.  To  provide  the  reader  with  an  accuracy  and  timing  study  of  FNOLl 
and  FN0L2  on  sample  problems. 

II.  Description  of  the  Method 

FN0L2  employs  the  Runge-Kutta  and  Adams-Moulton  fourth  order  techniques 
for  the  integration  of  a  system  of  first  order  ordinary  differential  equations, 
references  (1)  and  (2). 

A.  Form  of  Equations 

Let  the  system  of  equations  to  be  solved  be  given  in  the  form: 

(1)  y[  =  y^,  y2,  ...»  yjj)  i  =  1,2, ...,N 

yi<*0>  =  ^lo 


*In  most  cases,  adjustment  of  the  step  size  as  a  function  o'*  RELATIVE  error 
has  worked  very  well.  However,  in  cases  whero  at  least  one  of  the  Integrated 
variables  is  much  smaller  than  the  others  and  <  10-7  ,  this  method  has 
selected  an  unnecessarily  small  step  size.  In  this  event,  the  use  of 
absolute  truncation  error  adjustment  is  recommended. 
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Let  y.  be  the  value  of  y .  at  x  =  x  and  let  h  be  the  step  size  in 
in  1  n 

the  independent  variable  x. 

1.  Runge-Kutta  Formulas 

The  Runge-Kutta  method  uses  the  following  formulas  with  the 

appropriate  initial  conditions  to  go  from  the  n^^  to  the 
.st 

n+1  step. 


2. 


*  5  ^  i  '‘12> 

<*n  '>>  J'ln  ‘‘13’ 


J  ^‘‘12  ^13  ''14*  • 

Adame-Moulton  Formulas 

The  Adama-Moulton  predictor-corrector  formulas  for  the  system 
(1)  are: 


(2) 

V^p)  =  y  +  ]L 

^i,n+l  ^in  24 

(3) 

(c)  ,  h 

yi,n+l  =  ^in  ^  ^ 

,n-l  ^  ^i,n-2^ 


The  corrector  formula,  (3)  is  applied  only  once  so  that  only 
two  derivative  evaluations  are  needed  for  each  Adcuns-Moulton  integration  step. 
The  starting  values  needed  in  (2)  are  obtained  using the  Runge-Kutta  method 
(3  steps). 

Both  the  Runge-Kutta  and  the  Adame-Moulton  methods  incorporate 
round-off  control  features.  Tbxs  is  accomplished  hy  keeping  x  and  the  y^i^^ 
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in  double  precision  and  forming  the  sums  x  +  Zisx  and  ^  in  double 

precision.  All  FORTRAN  statements  are  evaluated  in  single  precision. 

B.  Automatic  Adjustment  of  Step  Size 

FN0L2  has  the  option  of  automatically  varying  the  step  size,  h,  to 
hold  the  Relative  or  Truncation  error  within  bounds  fixed  by  the 
user  (Adams-Moulton  Mode  only). 

The  Truncation  error,  T„  ,  is  approximated  by: 


the  Relative  error  is 


Procedure 

A.  Assignment  of  Variables 

In  order  to  integrate  N(1  <  N  <  30)  simultaneous  first  order 
ordinary  differential  equations  each  of  the  form: 

d/i 

j'i  =  dT  =  . ‘  .  " 

beginning  at  some  initial  point  and  ending  when 

. yN’J'i-yJ . J'n)  = 

one  proceeds  as  follows: 

For  notation  purposes,  assign  the  dependent  variables,  y^,  as 
elements  of  a  one-dimensional  array,  Y,  and  the  derivatives,  yj^. 


III. 

(4) 

(5) 
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as  elements  of  a  one-dimensional  array,  D.  Let 

1(1)  =  yi 

y(2)  =  72 


(6) 


Y(N)  =  y 
D(l)  =  y. 


D(N)  =  y‘ 

The  names  of  these  two  arrays  must  confora  to  the  usual 
FORTRAN  notation,  but  are  otherwise  arbitrary.  In  the  MAIN  program 
the  DIMENSION  of  the  Y  and  D  arrays  must  be  at  least  as  large  as 
(L  +  N  +  3).*  The  arrays  containing  the  dependent  variables  y^ 
and  its  derivatives  y|  should  not  be  placed  in  COMFCN  storage. 

One  or  more  additional  arrays  should  be  placed  in  COMMON  to  provide 
linkage  for  tables,  controls,  extraneous  calculations,  etc. 

6.  MAIN  Program 

The  MAIN  Program  must  assign  initial  values  to  the  dependent 
variables  (y^)  and  to  the  independent  variable  (x).  Somewhere 
in  the  MAIN  program  an  F  card  of  the  following  form  must  be 
inserted: 

F  DERIV,  TERM.  OUT,  TERM  1.  etc. 

- 7g.-..  — -  ^ - * - * - 


*L  and  N  are  control  variables  and  are  defined  on  p.6. 
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The  CALLING  sequence  is; 

CALL  FNOL2  (J,N,G,L,M,NE,X,y,D,DERIV, TERM, OUT) 

The  names  on  the  F  card  may  be  any  name  that  the  programmer  desires 

so  long  as  they  correspond  to  the  CALLING  sequence. 

The  controls  in  the  CALLING  sequence  for  FN0L2  are  as  follows; 

J  =  1  Use  Runge-Kutta  method  of  integrating  throughout 
(truncation  error  is  not  calculated). 

=  2  Use  Runge-Kutta  method  of  integrating  for  first  4  intervals, 
Adams-Moulton  (predictor  corrector)  until  termination  con¬ 
dition  in  Subroutine  TERM  is  reached  -  then  terminate  with 

4  iterations  using  Runge-Kutta  again.  Truncation  error  is 
automatically  calculated  in  Adams-Moulton  procedures. 

=  3  Use  Runge-Kutta  method  calculating  the  truncation  error  by 
repeating  every  second  calculation  with  a  doubled  interval 
of  integration  —  reference  (1).  The  step  size  is  not 
adjusted. 

N  =  Number  of  simultaneous  differential  equations  to  be  solved. 
They  are  set  up  in  the  DERIV  subroutine  and  may  be  any 
number  up  to  30. 

G  =  First  interval  of  integration. 

L  =  Number  of  D' s  greater  than  N  to  be  printed  out  when  the 
enclosed  OUTPUT  routine  is  followed  (see  p.  24).  If  the 
user's  output  routine  is  used,  L  may  be  ignored. 

M  =  Print  frequency - the  niunber  of  integration  cycles 

between  printouts.  If  M  =  0,  then  printing  is  determined 
by  values  assigned  to  Y(N  +  1)  and  y(N  +  2),  where  y(N  +  l) 
is  set  equal  to  some  running  variable  like  X,  D(l),  Y(l)  etc., 
and  Y(N  +  2)  is  a  constant  interval  in  y(N  +  1)  between 
printing  cycles.  This  feature  is  particularly  important  if 
the  variable  interval  of  integration  option  is  to  be  used 
since  printout  frequency  can  become  somewhat  meaningless 
when  the  interval  of  integration  changes  from  step  to  step. 

NE  =  Control  for  the  Interval  of  integration.  The  Interval  is 

increased,  let  alone,  or  decreased  depending  on  whether 

the  estimated  relative  or  truncation  error  is  less  than 

t„-NE-3  V.  ♦  irrNE-3  .  nn-NE 

10  ,  between  10  and  10  ,  or  greater  than 

10~”  .  Any  decrease  in  the  interval  of  integration  which 
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places  the  relative  or  truncation  error  below  the  band  is 
always  followed  by  an  advance  of  one  integration  step. 

Using  this  interval  before  a  subsequent  Increase  allows  the 
subroutine  to  proceed  through  difficult  regions  where  end¬ 
less  hunting  would  otherwise  result. 

NE  =  0  No  adjustment  of  G,  the  original  interval  of  integration 
takes  place,  and  Y(N  +  3)  is  ignored.  NE  must  be  zero 
when  J  =  1. 

y(N+3)*  =  -1.  And  IF  NE  ^  0  automatic  adjustment  of  the  interval  of 
integration  is  a  function  of  truncation  error. 

=  0.  And  IF  NE  ^  0  automatic  adjustment  of  the  Interval  of 

integration  is  a  function  of  the  Relative  truncation  error. 

X  =  Independent  variable.  The  initial  value  must  be  specified 
in  the  main  program  before  calling  FN0L2. 

7  =  The  name  given  to  the  array  which  contains  the  Dependent 
variables.  Initial  values  must  be  specified  in  the  main 
program  before  calling  FN0L2. 

D  s  The  name  given  to  the  array  which  contains  the  derivatives. 

The  first  N  D's  are  the  computed  values  of  the  differential 
equations  evaluated  at  current  values  of  X  and  7,  and  ob¬ 
tained  from  subroutine  DERIV.  The  D's  after  N  can  be  used 
to  print  out  auxiliary  quantities  or  compute  secondary 
qiiantitles  of  Interest,  except  in  subroutine  DERIV. 

DERIV  The  name  of  the  subroutine  in  which  the  differential  equations 
are  to  be  found.  In  multiple  phase  problems  this  name  can 
be  varied  as  DERIVl,  DERIV2,  etc.,  and  so  called  from  the 
main  routine. 

TERM  The  name  of  the  subroutine  in  which  the  termination  condi¬ 
tion  is  found.  Any  variable  or  combination  of  variables 
in  X,  7,  or  D  may  be  used  as  the  termination  condition. 

Several  termination  conditions  can  be  written  and  named 
TERMl,  TERM2,  etc.,  to  correspond  to  different  phases  of  a 
given  problem. 

OUT  A  separate  standard  output  routine  is  provided.  Any  output 
routine  desired  can  be  written  in  standard  FORTRAN  procedure 
however . 

The  names  of  the  3  subroutines  above  are  completely 
arbitrary.  The  only  restriction  is  that  they  must  corres¬ 
pond  to  the  F  cau^i  and  the  CALLING  sequence  for  FN0L2. 

"NOTE:  Y(N-»-3)  MUST  be  defined  in  the  executive  routine. 
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C .  Subroutine  DERIV 

Construct  a  subroutine  which  evaluates  the  derivatives,  yj^,  and 
places  the  results  in  the  D-array.  The  name  of  this  subroutine  is 
arbitrary.  The  first  statement  of  this  subroutine  must  be  of  the 
foi*a: 

SUBROUTINE  DERIV  (X,  Y,  D) 

No  recursive  formulas  should  be  evaluated  in  DERIV  since  it  may  be 
entered  a  variable  number  of  times  at  each  step,  depending  on  the 
mode  of  integration.  No  D  greater  than  N  should  be  used  in 
SUBROUTINE  DERIV. 


D.  Subroutine  TERM 

Construct  a  subroutine  which  evaluates  the  condition  and/or 
conditions  of  termination.  The  subroutine  must  be  of  the  form 
SUBROUTINE  TERM  (X,  Y,  D,  F) 

where  F  is  the  variable  which  determines  termination. F  may  be  equal 
to  f(X,  Y,  D).  Termination  occurs  when  F  =  0  or  undergoes  a  change 
of  sign.  Interpolation  and/or  iteration  is  automatic. 

Before  proceeding  to  a  discussion  of  the  output  routine  we 
digress  here  for  a  moment  to  give  an  example  of  the  above  two  sub¬ 
routines.  Suppose  that  it  is  desired  to  integrate  the  system 


beginning  at  some  Initial  point  (which  we  ignore  for  the  moment)  and 
terminating  when  yj^  =  0.  The  derivative  and  termination  subroutines 
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for  the  problem  are: 

SUBROUTINE  DERIV  (X,  Y,  D) 

DIMENSION  Y(l),  D(l) 

D(l)  =  Y(2) 

D(2)  =  -Y(l) 

RETURN 

SUBROUTINE  TERM  (X,  Y,  D,  F) 

DIMENSION  Y(l),  D(l) 

F  =  D(l) 

RETURN 

Note  that,  since  Y  and  0  are  arrays,  they  must  be  included  in  a 
DIMENSION  statement  in  each  subroutine. 

E.  Subroutine  OUTPUT 

A  standard  OUTPUT  subroutine  is  included  with  the  binary  deck  of 
FN0L2.  The  FORTflAN  listing  is  in  Appendix  A.  The  format  is  as 
follows: 


X  D(l) 

Y(l) 

truncation  error 

H( interval  of 

integration) 

D(2) 

Y(2) 

II 

D(3) 

D(4) 

D(5)  D(6) 

D(7) 

D(8) 

D(9) 

•  •  •  • 

.  .  D(L+N) 

Any  quantity  can  be  labeled  by  a  D  beyond  the  last  D  used  as  a 
differential  equation  (up  to  D(L+N))  and  it  will  be  printed  out 
automatically  in  the  above  format.  This  OUTPUT  routine  may  be 
replaced  by  any  FORTRAN  output  routine  with  the  following  CALLING 
sequence  and  DIMENSION  statement: 

SUBROUTINE  OUT  (X,  Y,  D,  ERROR,  N,  L,  H) 

DIMENSION  Y(l),  D(l),  ERROR  (1) 
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IV. 


Examples 

A.  Programming  Examples 
Example  No.  1 


d  V 

Use  the  Adams-Moulton  method  with  A  x  =  .01  to  integrate  — ^  = 

dx* 

and  terminate  when  x  =  11/2  . 

The  second  order  equation  must  be  reduced  to  two  first  order 
equations  of  the  following  form: 


-y 


dyi 

dF  -^2 

dy2 

dF  =  'yi 


(1) 

(2) 


The  initial  conditions  are; 

7^{0)  =  0 
y2(o)  =  1 

0  <  x  <  V2 

Using  the  enclosed  Output  routine  (Appendix  A)  write  the  following 
quantities  on  the  output  tape  at  every  fifth  step  of  integration: 

72*  yi*  yj’  ®(yi)»  '^72^  • 

Solution: 

For  FORTRAN  notation  we  make  the  following  correspondence: 

Y(l)  =  y^,  Y(2)  =  y2,  D(l)  =  y^,  D(2)  =  y^,  X  =  x. 

The  FORTRAN  program  is: 
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FOR 

C  MAIN  PROGRAM 

DIMENSION  Y<S| tOiS) 

X>0*0 
Y( 1)«0.0 
Y(2)>1.0 

F  DERIV. TERM, OUTPUT 

CALL  FNOL2  ( 2,2 »*T)1,0»5,0»X»Y,0, DERI  V»TERM, OUTPUT) 
STOP* 

END 


FOR 

SUBROUTINE  DERI V(X,Y.,0) 
DIMENSION  Y<ll ,D(1) 
D(1)>Y(2) 

D(2)—Y(l) 

RETURN 

END 


FOR 

SUBROUTINE  TERM(X,Y,0,F) 
DIMENSION  Y(1),0(1) 
F*X-1, 5707963 

return 

ENO 


*ST0P  slMuld  be  r«pl«e«d  vltb  aagr  laglUmt*  Slpataa  •ad.t* 
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Example  No.  2 

Integrate  the  equation  given  in  Example  1  with  the  Runge-Kutta 

method.  Read  In  all  controls  as  data.  Integrate  from  0  <  x  <  -^ 

with  h  =  .01  =  constat.  Terminate  at  it/2  and  use  the  Adams- 

Moulton  technique  with  a  variable  step  size  from  x/2  to  it  to  hold 

the  REUTIVE  Truncation  error  lO”®  <  RE^  <  10“^  (i.e.,  NE  =  5). 

2  2 

In  addition,  compute  +  y^  and  print  out  the  results  at  approxi¬ 
mately  equal  increments  of  0.1  in  y^. 

The  FORTRAN  program  is: 
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FOR 

C  MAIN  PROGRAM 

DIMENSION  Y(7) *D(7) 

1  READ  2»JfN.G.L*M»NE 

2  FORMAT! 213, F3. 3. 313) 

3  IF(J)  4,4,5 

4  STOP 

5  X=0.0 
Y(1)=0.0 
Y(2 ) =1.0 
Y(4) =0.1 

F  DERIV.TERMl ,T£RM2, OUTPUT 

CALL  FNOL2( J»N,G,L, M,NE,X,Y,D, DERI V,TERM1, OUTPUT) 

READ  2,NE',J 

Y(N+3)=0.0 

CALL  FNOL2(  J»N,G»  L,M,NE,X»Y,  D;»  DERI  V,TERM2,  OUTPUT) 
GO  TO  1 
END 
FOR 

■SUBROUTINE  DERIV(X,Y,D) 

DIMENSION  Y(1),0(1) 

D  (  1 )  =  Y  {  2  ) 

D{2)=-Y(l) 

Y(3)=Y(1) 

RETURN 

END 

FOR 

SUBROUTINE  TERMl ( X , Y ,D, F ) 

DIMENSION  Y(l) ,D(1) 

F  =  D( 1 ) 

D{3)=0( 1 )**2+D( 2)**2 
RETURN 
END 
FOR 

SUBROUTINE  TERM2 ( X , Y, D , F ) 

DIMENSION  Y ( 1 ) ,D< 1 ) 

F=X-3. 1415927 

D(  3) =D( 1 )**2+D( 2)**2 

RETURN 

END 


THE  DATA  FOR  THE  PROBLEM  ARE 


CARD 


j  N  G  L  M  NE 
3  2.011  0  0 


NE  J 

CARD  2  5  2 


CARD  3  0 
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Example  No.  3 

Compute  the  altitude  and  velocity  of  a  point  mass  dropped  vertically 
from  an  altitude  H.  Assume  that  IBAG  and  GRAVITY  are  the  only 
forces  acting  on  the  body.  Read  in  all  controls  and  initial  con¬ 
ditions  as  data.  Print  out  VELOCITY,  ALTITUDE  and  MACH  NO.  EVERY 


10  TIME  STEPS.  Use  an  initial  step  size  of  0.1  sec.  and  hold 

Terminate  when  H  <  H  TERM. 

The  equation  of  motion  is: 


lO”*^  <  RE^  <  10“^ 


±-Z  = 

dt^ 


_  -dV 
=  dt  = 


Cn  P  A 

D  ^o 
2M 


e-py  v2 


where 

2 

g  =  32.174  =  constant  =  acceleration  of  gravity  (ft/sec  ) 
Cjj=  .1  =  constant  =  drag  coeff.  dimensionless 

=  .0034  =  constant  =  density  of  air  Ib/fv 
P  =  22^5  =  constant 

M  =  Mass  =  constant  =  mass  of  body  (slugs) 

V  =  =  velocity  (ft/sec  ) 

y  =  H  =  altitude  (ft) 

A  =  1.0  =  area  of  body  (ft^) 

C  =  sound  speed  =  1116.4  (ft/sec)  constant 

H.  =  Termination  Alt.  (ft) 
term 


Let  y  =  P. 

•  •• 

p  =  y  =  -g 

and  y  =  P 


Then 
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become  our  two  first  oi*der  equations. 

Use  the  following  notation  for  programming  pvtrposes. 

Derivatives  Integrated  variables 

D(l)  =  y  C(l)  =  y 

D(2)  =  y  C(2)  =  y 


Data  and  Intermediate  Storage 


Y(l) 

=  V 

y(5)  = 

A 

Y(2) 

=  ^o 

Y(6)  = 

e-Py 

Y(3) 

=  Mass 

Y(7)  = 

_ (termination  Alt 

Y(4) 

Y(8)  = 

0  =  A 

Y(9)  = 

S  * 

Control  Variables 

K(l)  =  J  K(4)  =  M 
K(2)  =  N  K(5)  =  NE 
K(3)  =  L 


The  FORTRAN  program  is: 
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FOR 

C  MA*IN  PROGRAM 

OfMENSlON  Y{  10)  ♦0(5)  ,C(5) 

COMMON  YtK 
111  FORMAT{5I5) 

10  FORMAK 5E14,5) 

11  READ  lOf (Y( I )  il=lf8) 

IF(Y(2) ) 20^20  ♦  30 

F  DERIVfTERM^OUTPUT 

30  READ  111> (K( I ) ♦ Ial,b) 

C(1)*Y{1) 

C(2)=Y(2 ) 

N-K(2) 

C(N+3)=0.0 
T  =  0.0 

Y(9)  =  (Y(4)*Y(5)*.00  34)/2.*Y(3  ) 

CALL  Fi\0L2(K(l  )  ♦l<;(2)tY(8)^K(3)^K.(4)tK(5)  ♦  T  ♦  C  ♦  D  ♦  DER I  V  ♦  T  ERM  tOU  I  PUT  ) 
GO  TO  11 
20  STOP 
END 


FOR 

SUBROUTINE  DERIV(TfCtD) 
DIMENSION  Y(10)  ♦D(5)tC(5) tKIb) 
COMMON  YiK 

Y(6)*£XPF(-C(2) /22000. ) 

D( 1 ) =-32*174+Y( 9)*Y{6)*C( 1)**2 
D(2)=C(1) 

RETURN' 

END 


FOR 

SUBROUTINE  TERM ( T fC tO ♦ F » 
DIMENSION  Y(.10)  tDlS)  fCCS)  ♦IC(5) 
COMMON  YfK 
FaC{2)-rY(7) 

D(3)=C(1)/1116.4 

RETURN 

END 
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The  data  for  the 

problem  are; 

V 

yo 

Mass 

A 

card  1 

0.0 

100000. 

.7 

.1 

1. 

blank 

Term 

G 

card  2 

- 

10000. 

.1 

J 

N 

L 

M 

NE 

card  3 

2 

2 

1 

10 

4 

Results  of 

Test  Problems* 

1.  The  equations  in  Example  1 

dyi 

IT  - 

^2 

(1) 

and 

dx 

-^1 

(2) 

have  th«  solutions 

=  sine  x 
and 

y^  =  cos  X 

Equations  (1)  and  (2)  were  integrated  nunerically  from  0  <  x  <  3k 
using  FNOLl  and  FN0L2  with  a  constant  step  size  of  .01.  The  results 
of  the  Integration  and  an  error  analysis  are  presented  in  Tables  1 
and  2.  The  FNOLl  (sing}.e  precision)  solution  required  0.8  minutes 
machine  time  on  the  IBM  7090.  The  FN0L2  (double  precision) 
integration  took  0.9  minutes. 

^Unless  otherwise  stated,  all  integrations  were  performed  in  the  Adame- 
Moulton  mode. 
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2.  The  following  equation  which  Is  given  In  Hildebrand,  reference 
(3)  as  an  example  of  Instability  with  Milne's  method,  Is  In¬ 
tegrated  numerically  from  0  <  x  <  9,  with  constant  step  size 
of  .01. 

The  equation  Is: 

=  2  -  2y  ;  y(0)  =  2 

The  results  are  compared  to  the  true  solution  y(x)  =  e  +1  In 

Table  3  and  as  can  be  seen  the  numerical  solution  converges 

monotonlcally  to  the  true  solution.  The  Integration  required 

900  steps  and  .58  minutes  running  time  on  the  IBM  7090.  The 

—8 

maximum  error  was  1.49  x  10 

The  run  was  repeated  with  N£  =  4  and  It  required  only  51  steps 
and  .03  mins,  of  machine  time.  The  maximum  error  was  6.69x10""^. 
3<  Y  =  TAN  X  was  differentiated  twice  and  appropriate  substitutions 
were  made  to  arrive  at  the  following  second  order  equation: 
y  =  2  y  y 

The  above  equation  was  integrated  with  a  step  size  of  .01  frop 

0  <  X  <  1  with  y(0)  =  0  and  y(0)  =  1.  The  maximum  error  observed 
—8 

was  1.  X  10  .  The  rxm  required  0.6  min.  on  the  7090. 

The  run  was  repeated  with  NE  =  4  and  the  maximum  error  was 

3 >3  X  10^^  .  The  running  tine  was  reduced  to  .07  min. 

4«  a  =  uA  cos(wT)  was  integrated  from  0  <  T  <  -^  using  a  constant 

step  size  of  .1  with 

w  =  1 

A  =  1 
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The  maximum  absolute  error  (e^)  observed  was  2.9  x  lO”^ 

=  lAsinuT  -  /  a  | 

Using  the  same  At  and 
u  =  8 

A  =  15 

the  integration  becomes  unstable. 

The  above  equation  was  integrated  from  0  to  3n:  with  variable 
step  size  using  both  FNOLl  and  FN0L2.  The  .results  are  presented 
in  Table  4>  but  may  be  summarized  as  follows: 

The  FNOLl  (single  precision)  run  with  NE  =  6  required  1.58  min. 
machine  time  with  a  maximum  absolute  error  of  3.95  *  10  . 

The  FN0L2  (double  precision)  run  with  NE  =  4  took  1.0  min.  and 
had  a  maximum  error  of  4.13  x  lO”^. 


V.  Summary 

The  primary  objectives  of  this  report  have  been  to  present  the  user  with 
a  detailed  description  of  the  use  of  FN0L2  and  also  to  provide  a  comparison 
of  timing  and  accuracy  between  FNOLl  and  FN0L2. 

The  test  results  show  that  the  use  of  double  precision  in  FN0L2  does  not 
necessarily  increase  numing  time  when  a  variable  step  size  is  used.  Table  4 
shows  a  case  where  the  use  of  double  precision  has  produced  higher  accioracy 
with  reduced  running  time.  This  example  is  not  Intended  to  be  conclusive, 
since  an  optimal  step  size  for  a  given  problem  may  produce  equivalent  accuracy 
using  only  single  precision. 

The  automatic  adjustment  of  step  size  is  a  veiT-  valuable  feature.  It 
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can  result  in  significant  savings  of  nachlne  time  with  a  possible  reduction 
In  accuracy.  The  trade  of  accuracy  for  reduced  running  time  Is  a  choice 
which  the  user  must  make.  The  choice  should  be  based  on  careful  analysis  of 
the  problem,  experience  from  previous  results,  and  accuracy  requirements  for 
a  particular  problem. 
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iffBDn  A 

FOR 

SUBROUTINE  FN0L2 ( J .N ,G . L .M . NE » X . Y . D . DER I V , T ERM .OUTPUT ) 

C  6-63 

DIMENSION  Y( 50) .Dt 60) ,Yb(30.6 ) .GI2 ( 30 ) .GI 3 (30 )  .GI4(  30).EF(30). 
1  EFK  30  )  .EF2 (30).EF3(30).Y1(30) .ERROR (30) . HA ( 3 1 ) . Y A ( 50 ) . DA ( 50 ) » 
2YC ( 30) . YP( 30) . YAL ( 50 ) . YL ( 50 ) . YCL ( 30 ) 

ERASE  ( ( YB( I . J) .J=l .5 ) .ERROR( 1 ) . YL ( I ) , I = 1 . N ) . XU . NA . NF . NG . F , FA . 
iFB.FC.FO 
6  NB  =  1 

EC=Y(N+3) 

IF(EC)  1000.1.1 

1000  IF(EC+1. ) 1001.1,1001 

1001  EC=0. 

1  H  =  G 

2  HZ  =  H 

3  LN=N+XMAX0F(L,3) 

13  ENE=NE 

14  IF(J-2)  16,15,21 

15  IF(NE) 18,16,18 

16  JA=4 

17  GO  TO  22 

18  RE1=10.**(-ENE) 

19  RE2=10.#*(-ENE-3.0) 

20  RE(M*10.**(-EN,E-1.5) 

21  JA=1 

22  CALL  DERIV(X.Y,D) 

DO  300  1=1, N 
GI2( I )=D( I ) 

GI 3( I )=D(  I  ) 

GI4( I )=D(  I  ) 

300  EF( I )«D( I  ) 

27  CALL  OUTPUT(X,Y,D, ERROR, N,L.,H) 

28  FD=Y(N+1) 

29  IF(J-2)30,129,30 

30  GO  TO(31 ,37,35,37 ) ,JA 

31  DO  33  1  =  1, LN 

32  YA( I ) =Y( I ) 

YAL(I)  =YL(I) 

33  DA( I )=D( 1) 

34  GO  TO  37 

35  HB=H 

36  H=2.*H 

37  HD2=.5»H 

00  39  1=1, N 

38  YB( I ,NB)=0( I ) 

XL  =  D( I )  *  H02 

39  CALL  DPA( Y( I ) ,YL( I ) ,XL,0.,Y1( I ) ,XXL) 

CALL  DPA  (  X,XU,HD2,0,.  ,XX,XXL) 

40  CALL  DER-IV  ■  (XX,Y1  ,GI2) 

41  DO  42  1=1, N 

XL  =  GI2( I )*HD2 

42  CALL  DPA{ Y( I ) ,YL( I ) ,XL,0.»Y1( I ) ,XXL) 

CALL  DPA  (X,XU,H02,0,,XX,XXL) 

43  CALL  DERIV  (XX, Yl,  GI3) 

44  DO  45  1=1, N 
XL=GI3( I )*H 

45  CALL  DPA(Y( I ) ,YL( I ) ,XL,0.»Y1( I ) ,XXL) 

•a. 
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CALL  DPA( X,XU,H,0. .XX»XXL) 

46  CALL  DERI V( XX.Yl ,GI4) 

47  HD6  =H/6, 

GO  TO ( 48  * 5 5  * 60  *  66 ) * JA 

48  DO  52  1=1, N 

XL=(D(I)  +  2.#(GI2(I)  +  613(1))  +GI4(I))#HD6 

49  CALL  DPA  ( Y (  I  )  , YL ( I  )  ,XL  ,0. » Y ( I ) , YL ( 1 ) ) 

50  YC(  I.)=Y(  I  ) 

YCL(  I  ) =YL( I ) 

51  Yd  )  =YA(  I  ) 

YL(I)  =YAL(I) 

52  ERROR (I ) =0. 

53  JA=3 

54  GO  TO  35 

55  DO  57  1=1, N 

,-;l=(D(I)  +  2.»(GI2(I)  +  GI3(I))  +GI4(I))#HD6 

56  CALL  DPA  ( Y (  I  )  , YL (I  ) ,  XL  ,  0 . , Y ( I ) , YL ( I ) ) 

57  ERROR( I )  =  (Y( I )-YP( I  )  )/15. 

58  JA=1 

59  GO  TO  69 

60  DO  62  1=1, N 

61  Y (  I ) =YC (  I  ) 

YL(I)  =YCL(I) 

XL=(D(n  +  2.»(GI2(I)  +  GI3(1))  +GI4tI))#HD6 

62  CALL  DPA( YA { I )  ,YAL(  I  )  ,XL,0. ,YP( I )  ,XXL) 

63  H=HB 

64  JA=2 

65  GO  TO  69 

66  DO  68  1=1, N 

XL=(0(I)  +  2.#(GI2(I)  +  GI3(1))  +G14(I))#HD6 

67  CALL  DPA  ( Y (I  )  , YL (  I  ) , XL  ,  0 . , Y ( I )  ,  YL ( I ) ) 

68  ERROR! I ) =0. 

69  CALL  DPA  (  X  ,  XU , H , 0 . ,X , XU ) 

70  CALL  0ERIV(X,Y,D) 

71  FC=F 

72  CALL  TERM(X,Y,0,F  ) 

73-  IF(ABSF(F)-1.0E-5  )731,731,733 

731  NF=5 

732  GO  TO  124 

733  IF(F)74,124,76 

74  FA=1. 

75  GO  TO  77 

76  FB=1. 

77  IF(FA-FB)83.,78,83 

78  NF=NF+1 

79  JA=4 

80  NB=1 

81  H=H»F/(FC-F) 

82  IF(NF-4)37,37,124 

83  IF('NE)84,117,84 

84  IF(JA-l)  117,85,117 

85  IF(J-3)  86,117,86 

86  DO  95  1=1, N 

IF( Y( I ) )886,885,886 

885  HA ( I ) =1000. 

GO  TO  95 

886  IF(EC)  880,890*87 

87  IF( ABSF(Y( I ) )-EC)  880,880,890 

880  IF(ABSF(ERROR( I ) )-RE2)  882,94,881 

881  IF (ABSF (ERROR! I ) ) -REl ) 94 ,94 , 882 
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8  82  HA(  I  )  =H^^(  REM/ (AbSF(  ERROR  (  I  )  » .  OOOOOOOOOl )  )**(  .2) 

883  GO  TO  95 

890  IF(ABSF( ERROR ( I )/Y( I) > -RE2 ) 892 , 94 , 89 1 

891  IF(ABSF(ERROR( I l/Y( 1)  ) -RE  1 ) 94 , 94 , 892 

892  HA  (  I  )  =H*  (  REM/  (  ABSFI  ERROR (  D/Y  (  I  )  )■+. OOOOOOOOOl  ))'**(  .2  ) 

893  GO  TO  95 

94  ha ( I ) =H 

95  CONTINUE 

96  HA(N+1)=HA(N) 

HB  =  HA( 1 ) 

97  DO  98  I=2*N 

98  HB  =  MIN1F(HA( I ) ,HB) 

99  IF(H-HB) 100*117.101 

100  IF(HZ-H)101.101.116 

101  DO  103  1  =  1  ,LN 

102  Y ( 1 ) =YA{ I ) 

YLU  )  =YAL  (  I  ) 

103  D( I )=DA( I ) 

104  IF(NB-6) 107.105.105 

105  CALL  DPS  (X.XU.H,0..X,XU) 

106  GO  TO  109 

107  CALL  DPS  (X.XU,2.»H,0..X,XU) 

108  HZ=H 

109  H=HB 

CALL  DERIV(X.Y.D) 

110  NB=1 

111  XABS=ABSF{ ,0000ai*X) 

112  IF(ABSF(H)-XABS)113.113.117 

113  NG=NG+1 

114  H=SIGNF( XABS.HB) 

115  IFlNG-10) 124,126.126 

116  HZ=H 

117  IF(M) 118.118.121 

118  IF ( ABSF( ( Y( N+1 ) -FD) )-Y( N+2  >>29.119. 11° 

119  FD=Y(N+1> 

120  GO  TO  124 

121  NA=NA+1 

122  rF(M-NA) 123.123.29 

123  NA=0 

124  CALL  OUTPUT(X,Y.D, ERROR. N,L,H> 

125  IF(NF-4>29,29.126 

126  PRINT  127 

127  FORMAT(1HO> 

128  RETURN 

129  NB=NB+1 

130  IF(NB-6>30.131.136 

131  DO  134  1  =  1, N 

132  EF3( I >=YB( I  ,3> 

133  EF2( I >=YB( I ,4> 

134  EF.K  I  >=YB(  I  ,5> 

135  GO  TO  137 

136  NB=10 

137  HD24  =H/24. 

DO  138  I=1,N 

XL  =(55.*D(I>  -59.*EF1 ( I ).  +37.*EF2( I >  -9.»EF3 ( I > > *HD24 

138  CALL  DPA(  Y(  I  >  .YLU  >  .XL.O.  .YP(  I  >  »XXL> 

CALL  DPA  (X,XU.H.0..XX,XXL> 

139  CALL  DERIV(XX,YP.EF> 

140  DO  142  1=1. LN 

141  YA(I>=Y(I> 
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YAK  I)  =YL(I) 

142  DA{ I ) =D(  I  ) 

143  DO  148  1  =  1. N 

XL  =  (9.#EF(I)  +19.»D(I)  -5.*EF1(I)  +EF2 (  I  ) ) *HD24 

144  CALL  DPA(Y(  I  )  .'YL(  I  )  .XL.O..Y<  I  )  »yL(  I  )  ) 

145  €RROR( I )=( YP( I )-Y ( I ) )/14. 

146  EF3( I )=EF2( I ) 

147  EF2{ I )=EF1 ( I ) 

148  EFl ( I )=D(  I  ) 

149  GO  TO  69 

END 


FOR 

SUBROUTINE  OUTPUT ( X.Y.D. ERROR .N.L .H ) 
DIMENSION  Y(50) .D(50) »ERROR(30) 

PRINT  1 ,X,D( 1 1 .Y ( 1 » .ERROR! 1 ) »H 

1  FORMAT! 1H0.1P5E14. 7) 

IF!N-1)  20.20.10 

10  PRINT  2.  !  !D!  I  )  .Y!  I ) .ERROR! I ) ) . I=2.N) 

2  FORMAT ! 15X.1P3E14. 7) 

20  IF!L)50.50.30 

30  NZ=N+1 
NL=N+L 

PRINT  3.!D! I ) .I=NZ.NL) 

3  FORMAT! 15X.1P5E14, 7) 

50  RETURN 

END 
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TABLE  1 


F.NOL  1 

FNOL2 

TABLE 

S.P. 

D.P. 

X 

COSCX.) 

COS(X) 

COS(X) 

ERROR 

ERROR 

0.0 

1.0000000 

1.0000000 

1.0000000 

0. 

0.00000000 

0*5 

0.87758243 

0.87758256 

0.87758256 

0.00000013 

0.00000000 

1*0 

0.54030212 

0.54030234 

0.54030231 

0.00000019 

0.00000003 

1.5 

0.070737214 

.070737256 

0.070737200 

.00000012 

.00000005 

2.0 

-0.41614642 

-0.41614677 

-0.41614684 

.00000042 

.00000007 

2.5 

-0.80114267 

-0.80114356 

-0.80114362 

.00000095 

.00000006 

3.0 

-'0.98999119 

-0.98999248 

-0.98999250 

.00000131 

.00000002 

3.5 

-0.93645532 

-0.93645673 

-0.93645669 

.00000137 

.00000004 

4.0 

-0.65364260 

-0.65364373 

-0.65364362 

.00000102 

.00000011 

4.5 

-0.21079557 

-0.21079597 

-0.21079580 

.0000002-3 

.00000017 

5.0 

0.28366137 

-0.28366200 

0.28366218 

.00000081 

.00000018 

5.5 

0.70866791 

0.70866962 

0.70866977 

.00000186 

.OOOOOD15 

6.0 

0.96016772 

0.96017022 

0.96017029 

.00000257 

.00000007 

6.5 

0.97658492 

0.97658766 

0.97658763 

.00000271 

.00000005 

7.0 

0.75390015 

0.75390243 

0.75390225 

.00000210 

.00000018 

7.5 

0.34663444 

0.34663559 

0.34663532 

.00000088 

.00000027 

8.0 

-0.14549924 

-0.14549973 

-0.14550003 

.00000079 

.00000030 

8.5 

-0.60200945 

-0.60201164 

-0.60201190 

.00000145 

.00000026 

9.0 

-0.91112657 

-0.91113012 

-0.91113026 

.00000369 

.00000014 
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TABLE  2 


X 

FNOLl 

SIN(X) 

FN0L2 

SIN(X) 

TABLE 

SIN(X) 

S.F* 

ERROR 

o.p* 

ERROR 

0.0 

•  5 

0. 

0.47942,542 

0. 

0.47942552 

0. 

0.47942554 

0.000000000 

0.00000012 

0. 

0.00000002 

1.0 

0.84147059 

0.84147097 

0.84147098 

0.00000039 

*00000001 

l.S 

0.99749435 

0.99749498 

0.99749499 

0.00000064 

•00000001 

2.0 

0.90929670 

0.90929745 

0.90929743 

0.00000073 

•00000002 

2.3 

0.59847157 

0.59847221 

0.59847214 

0.00000057 

•00000007 

3.0 

0.14111992 

0.14112012 

CL.14112001 

0.00000009 

•00000011 

3.5 

-0.35078259 

-0.35078310 

-0*35078323 

.00000064 

•00000013 

4.0 

-0.75680106 

-0.75680239 

-0.75680250 

.00000144 

.00000011 

4.5 

-0.97752814 

-0.97753007 

-0.97753012 

.00000198 

.00000005 

5.0 

-0.95892226 

-0.95892433 

-0.95892428 

.00000202 

.00000005 

5.5 

-0.70553880 

-0.70554047 

-0.70554033 

.00000153 

•00000014 

6.0 

-0.27941496 

-0.27941572 

-0.27941550 

.00000054 

•00000022 

6.5 

0.21511917 

0.21511974 

0.21511999 

.00000082 

•00000025 

7.0 

0.65698439 

0.65698639 

0.65698660 

.00000221 

•00000021 

7.5 

0.93799679 

0.93799988 

0.93799998 

.00000319 

.00000010 

e.o 

0.98935486 

0.98935829 

0.98935825 

.00000339 

.00000004 

e.5 

0.79848440 

0.79848731 

0.79848711 

.00000271 

.00000020 

9.0 

0.41211718 

0.41211880 

0.41211848 

.00000130 

•00000032 

26- 


NOLTR  63-171 


TABLE  3 

y  =  /  i  =  e”^  +  1  ;  y(0)  =  2 


X 

Y(x)^ 

|YC  -  1 

0.0 

.5 

2.0 

1.3678795 

1.3678794 

6.00  X  10"® 

1.0 

1.1353353 

1.13533528 

1.2  X  10"® 

1.5 

1.0497871 

1.04978707 

3  X  10"® 

2.0 

1.0183156 

1.01831564 

6  X  10"® 

2.5 

1.0067379 

1.00673794 

4  X  10"® 

3.0 

1.0024787 

1.00247875 

5  X  10"® 

3.5 

1.0009119 

1.00091188 

1.2  X  10"® 

4.0 

1.0003355 

1.00033546 

4  X  10"® 

4.5 

1.0001234 

1.00012339 

1  X  10"® 

5.0 

• 

1.0000454 

• 

1.0000455 

1  X  10”® 

• 

• 

8.39 

9.0 

• 

• 

1.00000000 

1.0000000 

1.00000000 

1.00000000 

*Table  Values 
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TABLE  4 

e  =  (Asin(uT)  -  /  A 

ucos(ut) 1 

t 

FNOL  1 

FN0L2 

e 

t 

0.0 

0.0 

0.0 

1.0 

3.48  X  10“^ 

7.08  X  10“^ 

2.0 

1.12  X  10“^ 

1.70  X  lO”^ 

3.0 

3.94  X  lO”^ 

2.17  X  10'^ 

4.0 

7.05  X  10“^ 

3.0  X  10"^ 

5.0 

1.25  X  lO”^ 

3.83  X  10"^ 

6.0 

1.76  X  lO"^ 

4.13  X  10”^ 

7.0 

2.33  X  10“^ 

3.80  X  10“*^ 

8.0 

3.24  X  10”^ 

3.33  X  10“^ 

9.0 

3.95  X  10”^ 

1.43  X  10“^ 

3  nr 

3. a  X  10"^ 

3.06  X  10”^ 

Machine  Tlote 

1.58  Min. 

1.0  Min. 
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